function hfile = voi_AddSphericalVOI(hfile, c, r)
% VOI::AddSphericalVOI  - add a spherically shaped VOI to the object
%
% FORMAT:       [voi = ] voi.AddSphericalVOI(c, r);
%
% Input fields:
%
%       c           1x3 coordinate (center)
%       r           1x1 radius (must be > 0 and < 128)
%
% Output fields:
%
%       voi         VOI with added VOI (integer coordinates only!)

% Version:  v0.7g
% Build:    9040105
% Date:     Apr-01 2009, 5:38 AM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin ~= 3 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'voi') || ...
   ~isa(c, 'double') || ...
    numel(c) ~= 3 || ...
    any(isinf(c) | isnan(c) | c < -128 | c > 256) || ...
   ~isa(r, 'double') || ...
    numel(r) ~= 1 || ...
    isinf(r) || ...
    isnan(r) || ...
    r < 0 || ...
    r > 128
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);

% coordinate valid in current system
if (strcmpi(bc.ReferenceSpace, 'tal') && ...
    any(c > 128)) || ...
   (~strcmpi(bc.ReferenceSpace, 'tal') && ...
    any(c < 0))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid coordinate given (wrong coordinate system).' ...
    );
end
if strcmpi(bc.ReferenceSpace, 'tal')
    bb = [-128, 128];
else
    bb = [0, 256];
end

% create grid
rc = round(c);
cr = ceil(r + 0.5);
[xg, yg, zg] = ndgrid( ...
    rc(1)-cr:rc(1)+cr, rc(2)-cr:rc(2)+cr, rc(3)-cr:rc(3)+cr);
xg = [xg(:), yg(:), zg(:)];
xg(any(xg < bb(1), 2) | any(xg > bb(2), 2), :) = [];
yg = sqrt(sum((xg - c(ones(1, size(xg, 1)), :)) .^ 2, 2));
xg(yg > r, :) = [];

% create VOI structure
voi = struct( ...
    'Name', sprintf('Sphere_%.1f_%.1f_%.1f_-_r%.1f', c(1), c(2), c(3), r), ...
    'Color', floor(255.99 * rand(1, 3)), ...
    'NrOfVoxels', size(xg, 1), ...
    'Voxels', xg);

% set back
bc.VOI(end + 1) = voi;
bc.NrOfVOIs = numel(bc.VOI);
bvqxfile_setcont(hfile.L, bc);
